% normals = getNormals(lights,grayImgs,grayMask)
%
% input:
%   * lights - A 3-by-numImages matrix representing the lighting direction
%   of each image in the image series. lights(:,i) is the lighting
%   direction vector of the ith image in the series. A lighting direction
%   vector is of the form (l_x, l_y, l_z), where the x axis is along the
%   column index increasing direction (towards right) and the y axis is
%   along the row index increasing direction (downwards). The z axis points
%   out of the image.
%   * grayImgs - A numRows-by-numColumns-by-numImages matrix representing a
%   series of images taken for the same scene under different lighting
%   conditions. grayImgs(:,:,i) is the ith image in the image series, 
%   presented in gray scale.
%   * grayMask - A numRows-by-numColumns matrix masking the scene related
%   pixels. However, this matrix doesn't consist of only 0 and 1 indicating
%   exactly whether a pixel is or is not a scene point. The values kind of
%   range from 0 to 255. In this method, if grayMask(i,j) is larger than
%   127.5, then pixel(i,j) in the image is considered to be a scene point;
%   otherwise, the background pixel.
%
% output:
%   * normals - A numRows-by-numColumns-by-3 matrix representing the
%   normals of each pixel in a numRows-by-numColumns image. normals(i,j,:)
%   is the normal vector of the scene point corresponding to pixel(i,j) in
%   the image. The normal vector is in the form of (n_x, n_y, n_z), where
%   the x axis is along the column index increasing direction (towards
%   right) and the y axis is along the row index increasing direction
%   (downwards). The z axis points out of the image.
%
% description:
% This method takes the lighting direction for each image in the image
% series, the image series (in gray scale), and the image scene point mask
% (in gray scale) and calculates the normal vector for each scene point.
% The normal vector is in the form of (n_x, n_y, n_z), where the x axis is
% along the column index increasing direction (towards right) and the y
% axis is along the row index increasing direction (downwards). The z axis
% points out of the image.

function normals = getNormals(lights,grayImgs,grayMask)
[nRs, nCs, nImgs] = size(grayImgs);
normals = zeros(nRs, nCs, 3);
threshold = max(max(grayMask))/2;
for i = 1:nRs
    for j = 1:nCs
        if grayMask(i,j) >= threshold
            intensities = double(reshape(grayImgs(i,j,:),1,nImgs));
            w = weight(intensities)';
            A = -repmat(w,1,3).*lights';
            B = w.*intensities';
            n = A\B;
            n = n/sqrt(n'*n);
            if sum(isnan(n)) ~= 0
                disp('Warning: wrong normal!');
            end
            normals(i, j, :) = reshape(n,1,1,3);            
        end
    end
end

end

function w = weight(intensity)
    w = 1-abs(127.5-intensity)/127.5;
    %w = intensity;
end

